Statistique et apprentissage pour la prévision de séries temporelles
Prévision des locations de vélibs à Londres
Diaby Drame & Cynthia Lepere
Janvier 2023
Introduction
Nous sommes deux étudiants habitués aux métros parisiens, surchargés tout autant que le trafic routier, comme dans d’autres grandes villes. Nous pouvons donc affimer que se déplacer en vélo présente beaucoup d’avantages…mais malheureusement, un inconvénient existe : parfois il n’y a plus de vélos disponibles dans une station désirée, ou alors plus de places pour les rendre.
Nous avons lu un article d’un enseignant chercheur en mécanique des fluides qui a créé justement un logiciel de prédiction, intégré dans l’application «La Bonne Station» de Keolis, destinée aux utilisateurs du VCub a Bordeaux (équivalent des velibs a Paris) et nous avons trouvé intéressant cet enjeux de prévision. Son objectif est de maximiser l’efficacité des opérations de régulation des vélos: Quelles stations approvisionner en priorité ? Quelles stations agrandir en priorité ? Où positionner de nouvelles stations ? Quelles stratégies adopter pour être le plus efficace possible ?
Présentation des données
Pour notre projet d’analyse prédictive reliée au vélibs, nous allons
travailler sur le jeu de données London bike.
Nous l’avons trouvé sur le site Kaggle.
Il correspond a des données de locations de vélos dans la ville de
Londres.
| timestamp | cnt | t1 | t2 | hum | wind | weather | holiday | weekend | season |
|---|---|---|---|---|---|---|---|---|---|
| 2015-01-04 00:00:00 | 182 | 3.0 | 2.0 | 93.0 | 6.0 | 3 | 0 | 1 | 3 |
| 2015-01-04 01:00:00 | 138 | 3.0 | 2.5 | 93.0 | 5.0 | 1 | 0 | 1 | 3 |
| 2015-01-04 02:00:00 | 134 | 2.5 | 2.5 | 96.5 | 0.0 | 1 | 0 | 1 | 3 |
| 2015-01-04 03:00:00 | 72 | 2.0 | 2.0 | 100.0 | 0.0 | 1 | 0 | 1 | 3 |
| 2015-01-04 04:00:00 | 47 | 2.0 | 0.0 | 93.0 | 6.5 | 1 | 0 | 1 | 3 |
Il contient \(17414\) observations, recensées sur deux ans, pendant chaque heure entre le 04/01/2015 à 00h00 et le 03/01/2017 à 23h00, et \(10\) variables :
- timestamp = l’horodatage
- cnt = le nombre de vélos empruntés
- t1 = la température réelle (en C°)
- t2 = la température ressentie (en C°)
- hum = l’humidité (en %)
- wind = la vitesse du vent (en km/h)
- weather = la météo (1=clair,2=peu nuageux,3=nuages fragmentés,4=nuageux/couvert…,7=pluie)
- holiday = jour de vacances(1) ou non vacances(0)
- weekend = jour de week-end(1) ou non week-end(0)
- season = la saison (0-printemps ; 1-été; 2-automne; 3-hiver)
\(\underline{Problématique}\)
Nous commencerons par chercher a prédire le nombre d’emprunts de vélos
des 4 derniers mois du jeu de données, avec des modèles qui utilisent
juste les variables explicatives.
Puis nous ferons un modele prédictif, afin de prevoir les futurs
emprunts de vélos.
Notre variable d’intérêt sera donc le nombre de vélos empruntés
«cnt».
Et nos variables explicatives seront les autres.
\(\underline{Traitement~des~données}\)
Nous commençons par verifier que nous n’avons pas de données
manquantes:
C’est bien le cas, donc nous n’avons pas besoin de faire d’imputation de celles-ci!
Ensuite, nous ajoutons plusieurs variables qui nous seront utiles pour la suite:
- day_of_month = le jour du mois
- month = le mois (1=janvier…12=decembre)
- day_of_week = le jour de la semaine (1=Dimanche…7=Lundi)
- hour = l’heure
- Posy = la position au cours de l’année
- date = la date sous la forme année-mois-jour
- Year = l’année
| timestamp | cnt | t1 | t2 | hum | wind | weather | holiday | weekend | season | day_of_month | month | day_of_week | hour | Posy | date | Year |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2015-01-04 00:00:00 | 182 | 3.0 | 2.0 | 93.0 | 6.0 | 3 | 0 | 1 | 3 | 4 | 1 | 1 | 0 | 0.0082418 | 2015-01-04 | 2015 |
| 2015-01-04 01:00:00 | 138 | 3.0 | 2.5 | 93.0 | 5.0 | 1 | 0 | 1 | 3 | 4 | 1 | 1 | 1 | 0.0082418 | 2015-01-04 | 2015 |
| 2015-01-04 02:00:00 | 134 | 2.5 | 2.5 | 96.5 | 0.0 | 1 | 0 | 1 | 3 | 4 | 1 | 1 | 2 | 0.0082418 | 2015-01-04 | 2015 |
| 2015-01-04 03:00:00 | 72 | 2.0 | 2.0 | 100.0 | 0.0 | 1 | 0 | 1 | 3 | 4 | 1 | 1 | 3 | 0.0082418 | 2015-01-04 | 2015 |
| 2015-01-04 04:00:00 | 47 | 2.0 | 0.0 | 93.0 | 6.5 | 1 | 0 | 1 | 3 | 4 | 1 | 1 | 4 | 0.0082418 | 2015-01-04 | 2015 |
Analyse descriptive
Voici le nombre de vélos empruntés au cours du temps:
Nos données ne recouvrent que deux ans, mais nous pouvons voir
notamment une saisonnalité.
Nous n’avons pas de changement signification d’une année à l’autre dans
la locations des vélos.
On remarque a priori qu’ils sont plus utilisés pendant l’été que
l’hiver.
Nous voyons deux valeurs anormales / deux pics élévés durant l’été 2015.
Les deux pics correspondent à des jours de grève des transports.
Regardons tout d’abord le nombre de vélos loués en fonction des différentes variables.
\(\underline{Variables~calendaires}\)
\(\bullet~Graphiques~annuel\)
Sur le premier graphique, nous voyons le nombre de locations de vélos à chaque mois tout au long de l’année. La courbe noire correspond a la premiere année de notre jeu de données (2015/2016), et celle en bleue a la deuxieme année (2016/2017).
Nous voyons un pic élévé au milieu des deux années (7eme mois: juillet) ce qui correspond aux étés. Tandis que les niveaux les plus bas sont en hiver (1er mois: janvier) et sont inférieur à ceux des autres saisons.
Cela nous confirme donc que les vélos sont moins empruntés pendant l’hiver.Sur le deuxieme graphique, nous avons la moyenne des deux années, et elle nous indique les même conclusions.
\(\bullet~Graphiques~hebdomadaire\)
Le premier graphique représente le nombre de locations de vélos par jour(1=Dimanche…7=Lundi) en 1 semaine.
Nous remarquons que les vélos sont plus loués pendant la semaine que pendant le week-end.Nous avons également représenté ces resultats sous la forme d’un diagramme à barres.
\(\bullet~Graphiques~quotidien\)
Le premier graphique montre le nombre de locations de vélos par heure pendant 1 jour.
La courbe bleue correspond aux jours de vacance, et celle en noire aux autres jour.
Le deuxieme graphique nous montre le nombre de locations de vélos par heure, pour chaque jours de la semaine.Nous voyons clairement (et heureusement!!) que les locations pendant les non vacances, correspondent aux locations durant les jours de semaines. Et que les locations pendant les vacances, sont similaires a celles des jours de week end.
Par rapport aux vacances : nous remarquons qu’ils sont plus loués aux alentours de 14heures.
Alors que le reste de l’année : ils sont plus loués aux moment des heures de pointe (8 et 18 heures)
Et de même par rapport aux jours : En semaine, le nombre total de locations est faible de 00h00 à 05h00 et augmente jusqu’au premier pic à 08h00, puis les déplacements se stabilisent et restent à peu près constants jusqu’à 15h00.
Le deuxième pic est atteint à 17h00 et continue ensuite à baisser jusqu’à environ 500 à 23h00.
Cela suggère potentiellement que les cyclistes utilisent les vélos publics pour se rendre au travail vers 08h00 puis en revenir vers 17h00 en semaine.
Le week-end, les locations commencent à augmenter lorsque le soleil se lève, et diminuent le soir.
\(\underline{Variables~liées~a~la~météo}\)
Nous remarquons que les vélos sont le plus loués lorsqu’il fait plus
chaud, et qu’il y a moins d’humidité.
La vitesse du vent n’a pas vraiment l’air d’influer sur les locations
par contre.
\(\underline{Corrélations}\)
Sur la matrice de corrélation, nous remarquons que notre variable d’interet cnt est le plus corrélée avec les variables t1, t2, hum et hour.
Sur le deuxieme “graphique” la proximité des points est déterminée à l’aide du clustering multidimensionnel.
Les variables temperature réelle t1 et temperature ressentie t2 apparaissent plus proche et sont jointes par un chemin plus foncé, elles sont donc les plus fortement corrélées. Nous utiliserons donc par la suite uniquement l’une des deux.
Les variables day_of_month, day_of_week, holiday et weekend forment comme un petit groupe, elles “donnent les mêmes explications”. Les autres variables explicatives ne sont pas tellement corrélées entre elles.
\(\underline{Autocorrélations}\)
L’autocorrélation (ACF) d’une série fait référence au fait que dans une série temporelle, la mesure d’un phénomène à un instant t peut être corrélée aux mesures précédentes (au temps t-1, t-2, t-3, etc) ou aux mesures suivantes (à t+1, t+2, t+3, …). L’autocorrélation mesure essentiellement la similitude entre les observations en fonction du décalage horaire entre elles. Une série autocorrélée est ainsi corrélée à elle-même, avec un décalage (lag) donné.
Sur notre graphique ACF, l’autocorrélation est la plus forte pour les décalages 1 et 24. (Et 0 mais ça ne signifie rien car l’autocorrélation de la série avec elle-même, sans décalage, est forcément de 1!)L’autocorrélation (PACF) partielle permet de mesurer l’autocorrélation d’un signal pour un décalage k “indépendamment” des autocorrélations pour les décalages inférieurs.
Sur notre graphique ACF, on a une forte autocorrélation pour un décalage de 1, et de 24… Or l’autocorrélation pour le décalage de 24 pourrait très bien être une conséquence de celle existant pour le décalage de 1, c’est ce pourquoi nous regardons ensuite le graphique PACF. La valeur au lag 1 est bien plus forte que pour les autres, ce qui nous laisse penser que l’autocorrélation observée au décalage 24 était un effet résiduel de l’autocorrélation pour un décalage de 1.
Si la série était stationnaire, nous verrions pratiquement toutes les lignes dans les intervalles de confiance a 95% représentés en verts. Or, nous voyons que chaque pic est hors de ces lignes.
Modélisation
Predictions de la fin du jeu de données
Notre objectif est de développer des modèles basés sur l’ensemble de
données d’entraînement et de faire des prédictions du nombre horaire de
location de vélos à l’aide de l’ensemble de données de test.
Nous allons donc estimer des modèles sur des données d’entrainement et
faire des prévisions sur des données de test.
Nous commençons par divisé notre jeu de données en une partie pour
entrainer nos modèles, et une partie pour les valider.
Nous gardons les 4 derniers mois pour nos validations.
Nous avons 14488 observations pour nos tests, et 2926 pour nos validations.
Nous allons appliqué plusieurs modèles afin de prédire le nombre
total de locations de vélos par heure, pour les quatres dernier
mois.
Nous les évaluerons grâce a differentes valeurs :
- La racine de l’erreur quadratique moyenne (RMSE) définie par: \(RMSE = \sqrt{MSE} = \sqrt{\frac{1}{n} \sum_{i=1}^n
(y_i–\hat{y}_i)^2}\)
Elle s’exprime dans la même unité que la variable à prédire.
Plus elle est élevée, moins le modèle est performant.
- L’erreur absolue moyenne en pourcentage (MAPE) définie par: \(MAPE\) = \(\frac{1}{n} \sum_{i=1}^n \displaystyle\left\lvert
\frac{e_i}{y_i} \right\rvert\) = \(\frac{1}{n} \sum_{i=1}^n \displaystyle\left\lvert
\frac{ (y_i–\hat{y}_i) }{ y_i } \right\rvert\)
Elle s’exprime en pourcentage et elle peut se traduire par le pourcentage moyen d’écart entre la valeur prédite et la réalité.
De même, plus elle est élevée, moins le modèle est performant.
- Le temps d’exécution
Remarque:
- Nous avons choisi d’utiliser la MAPE et non la MAE (erreur absolue
moyenne), pour deux raisons:
La première est qu’un pourcentage est plus facile à comprendre, et la
seconde est que les scores MAPE peuvent être comparés entre différents
modèles. Cependant, si vos valeurs réelles auraient été proches de 0,
nous aurions choisis la MAE!!
- La MSLE (erreur logarithmique quadratique moyenne) réduit l’importance
des erreurs pour de grandes valeurs réelles, comme le fait deja la MAPE.
Elle ne s’exprime pas dans une unité simple, ce qui la rend difficile à
interpréter, ce pourquoi nous ne la choisissons pas non plus!
\(\underline{GLM}\)
RMSE = 1005.917 MAPE = 1.536 utilisateur système écoulé
0 0 0 Ce modèle est très rapide, mais le résultat n’est pas terrible..
\(\underline{SARIMA}\)
Series: bikeCount
ARIMA(5,1,0)
Coefficients:
ar1 ar2 ar3 ar4 ar5
0.3228 -0.3496 -0.0672 -0.0203 -0.0606
s.e. 0.0076 0.0080 0.0084 0.0079 0.0076
sigma^2 = 416609: log likelihood = -137366.9
AIC=274745.8 AICc=274745.8 BIC=274792.4
RMSE = 1328.121 MAPE = 4.406 utilisateur système écoulé
0.01 0.00 0.01 Le modèle SARIMA est une extension du modèle ARIMA.
Nous otenons une RMSE d’environ 1328 vélos, ce qui est encore moins bon
que pour notre modele de regression. On constate que le modèle est très
bon sur les données d’entraînement et moins bon sur les données de test.
On peut donc soupçonner des cas de sur-apprentissage.
\(\underline{Arbres~de~décision}\)
RMSE = 534.295 MAPE = 0.788 utilisateur système écoulé
0.24 0.02 0.25 Sur le premier graphique, nous pouvons visualiser l’arbre obtenu.
Trois de nos variables hour, t1 et weekend ont été utilisées, et 14 nœuds terminaux ont été générés.
Sur notre ensemble d’entrainement, nous voyons que la variable d’heure hour est la première variable choisie pour la décision. Le nombre de vélos empruntés attendu, lorsque hour<7 est de 197, et 29% des échantillons tombent dans cette feuille. Lorsque hour est inférieure à 7heure, moins de vélos sont utilisés, et lorsque hour est >7h et <20h, la température commence à entrer en jeu.Sur le deuxieme graphique, nous pouvons voir l’erreur de validation croisée (relative) pour chaque sous-arbre, du plus petit au plus grand.
Dans “les coulisses”, rpart() applique automatiquement une validation croisée pour élaguer l’arbre. L’axe x en bas correspond a la complexité de l’arbre.
L’axe x en haut correspond au nombre de nœuds terminaux, donc la taille de l’arbre.
L’axe y a l’erreur de convergence.
Un bon choix de cp pour l’élagage est souvent la valeur la plus à gauche pour laquelle la moyenne se trouve sous la ligne horizontale. Nous voyons que l’endroit optimal pour tailler l’arbre correspond a 12 nœuds.Notre RMSE pour l’ensemble de test est d’environ 534 vélos!!
la MAPE vaut environ 0.79, et le temps d’execution est de 0.118 secondes.Remarque: Les arbres de décision sont des algorithmes où les hyperparamètres sont essentiels en ce qui concerne leur performance. Par exemple, un nombre de profondeurs maximales (maxdepth) trop élevées peut entraîner un surajustement ou une variance élevée, et si ce nombre est trop faible cela peut entraîner a l’inverse un manque d’ajustement ou un biais élevé. Il faut donc faire attention car l’objectif est d’avoir les meilleurs resultats, mais l’arbre pourrait être suradapté.
\(\underline{GAM}\)
RMSE = 603.294 MAPE = 1.056 utilisateur système écoulé
0.04 0.00 0.05 Plus rapide (seulement 0.014 secondes), mais la RMSE est d’environ 603 vélos, et la MAPE d’environ 1.06 : c’est un peu moins bien que l’arbre !
\(\underline{Randomn~forest}\)
RMSE = 557.427 MAPE = 0.835 utilisateur système écoulé
0.16 0.11 0.14 Notre RMSE pour l’ensemble de test est d’environ 549 vélos, ce qui est un tout petit peu plus que pour l’abre de décision, mais ils ont la même valeur de MAPE, environ 0.79, et nos forets aleatoires sont plus rapide (0.028 secondes)
\(\underline{XGBoost}\)
RMSE = 372.815 MAPE = 0.718 utilisateur système écoulé
0.04 0.00 0.02 La RMSE vaut environ 373 vélos et la MAPE environ 0.718!!! Nous obtenons les meilleurs resultats, l’extreme gradient boosting dépasse tout le monde! De plus le temps d’execution est seulement de 0.002 secondes.
\(\underline{Perceptron}(MLP)\)
Entrainons le perceptron multi-couches en faisant varient le nombre de neurones, le nombre de couches afin de trouver un modèle qui prédit mieux.
La fonction RELU \(f(x)=x^{+}=max(0,x)\) sera notre fonction d’activation. On utilisera le package Keras.
- cas d’une couche en entrée et sortie (pas de couche cachée).
128 neurones à la première couche et 1 à la dernière.
RMSE = 909.46 MAPE = 1.345 utilisateur système écoulé
0.44 0.82 0.32 - cas d’une couche en entrée, une couche cachée et une couche de sortie.
500 neurones à la première et deuxième couche et 1 à la dernière.